# Script to run QAP
library(sna)
library(doParallel)
library(doRNG)
library(biglm)


# Load the followers adjacency matrix
y <- readRDS("processed_data/retweets_adjacencyMatrix.Rds")
# Load the predicting matrix 
load("processed_data/QAP_predicting_matrices.RData")

diag(y) <- NA

y_v <- c(y)


Var.names <-c("State_Similarity",
              "Party_Similarity",
              "Chamber_Similarity",
              "Gender_Similarity",
              "Race_Similarity",
              "Difference_in_Legislatures_Profeshionalism",
              "Dem_Sender_Effect",
              "Rep_Sender_Effect",
              "House_Sender_Effect",
              "Female_Sender_Effect",
              "Profesh_Sender_Effect",
              "Black_Sender_Effect",
              "Latino_Sender_Effect",
              "Asian_Sender_Effect",
              "Mena_Sender_Effect",
              "Multi_Sender_Effect",
              "Native_Sender_Effect",
              "Democrat_Receiver_Effect",
              "Republican_Receiver_Effect",
              "House_Receiver_Effect",
              "Female_Receiver_Effect",
              "Profesh_Receiver_Effect",
              "Black_Receiver_Effect",
              "Latino_Receiver_Effect",
              "Asian_Receiver_Effect",
              "Mena_Receiver_Effect",
              "Multi_Receiver_Effect",
              "Native_Receiver_Effect",
              "Same_Party_Same_State",
              "Same_Chamber_Same_State",
              "Same_Gender_Same_State",
              "Same_Race_Same_State",
              "Contiguous_States")

for(i in 1:length(Var.names)){
  xi <- predicting_matrices[i,,]
  diag(xi) <- NA
  print(sum(is.na(xi)))
  if(i == 1){
    xdf <- cbind(c(xi))
  }
  if(i > 1){
    xdf <- cbind(xdf,cbind(c(xi)))
  }
}

send_id <- matrix(1:nrow(y),nrow(y),nrow(y))
rec_id <- t(send_id)
diag(send_id) <- NA
diag(rec_id) <- NA

xdf <- cbind(xdf,c(send_id),c(rec_id))

xdf <- data.frame(xdf)

names(xdf) <- c(Var.names,c("send_id","rec_id"))

xydf <- data.frame(y_v,xdf)

xynd <- na.omit(xydf)

covs <- paste(Var.names,collapse="+")

names(xynd)[1] <- "yij"

load("simulated_retweet_networks.RData")

for(j in 1:25){

  # qap unpermuted
  ys <- sim_amats[[j]]
  
  diag(ys) <- NA
  
  y_vs <- c(y)
  
  xydf$y_vs <- y_vs
  
  covs <- paste(Var.names,collapse="+")
  
  form <- as.formula(paste("y_vs~",covs,sep=""))
  
  system.time(est_true <-  biglm(form,data=xydf))
  
  # qap permuted
  
  set.seed(100416)
  perm_coef <- NULL
  
  niter <- 100
  for(i in 1:niter){
    y_p <- rmperm(ys)
    y_vp <- c(y_p)
    xydf$y_vp <- y_vp
    form <- as.formula(paste("y_vp~",covs,sep=""))
    print(system.time(est_perm <-  biglm(form,data=xydf)))
    gc()
    tstat <- summary(est_perm)$mat[,1]/summary(est_perm)$mat[,4]
    perm_coef <- rbind(perm_coef,tstat)
  }
  
  pgeq <- NULL
  for(i in 1:length(coef(est_true))){
    pgeq <- c(pgeq,mean(perm_coef[,i] > coef(est_true)[i]))
  }
  
  pleq <- NULL
  for(i in 1:length(coef(est_true))){
    pleq <- c(pgeq,mean(perm_coef[,i] < coef(est_true)[i]))
  } 
  
  pgabs <- NULL
  for(i in 1:length(coef(est_true))){
    pgabs <- c(pgabs,mean(abs(perm_coef[,i]) > abs(coef(est_true)[i])))
  } 
  
  save(list=c("perm_coef","pgeq","pleq","pgabs"),file=paste("./qap_power_results/qap_sim_retweet",j,".RData",sep=""))
  
}









